Pipeline

This was performed by merging only WT (n = 6), only TG (n = 6) and all the samples (n = 12) together at Isoseq3 Cluster, align to ERCC, cupcake collapse and chain three datasets, followed by SQANTI.

Pipeline for analysis

Classification file

# Read in sqanti classification file mapped to ERCC (post chaining WT and TG)
class <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/SQANTI/all_samples.chained_classification.filtered_lite_classification.txt", 
                    header = T, as.is = T, sep = "\t")

datatable(head(class), options = list(scrollX = TRUE))

Unique ERCC

# replace FL count with pb id, otherwise NA if 0
unique <- setNames(data.frame(matrix(ncol = 4, nrow = nrow(class))), c("isoform","FL.All_Merged","FL.TG_Merged","FL.WT_Merged"))
unique$isoform <- class$isoform
for (j in c("FL.All_Merged","FL.TG_Merged","FL.WT_Merged")){
  for(i in 1:nrow(class)){unique[[j]][i] <- ifelse(class[[j]][i] > 0, paste(class$chrom[i]),paste("NA"))}
}


# count the unique pbid (/ERCC) per dataset
p1 <- apply(unique[,-1], 2, function(x) length(unique(x))) %>% melt() %>% 
  rownames_to_column(., "dataset") %>%
  mutate(type = c("All","NotAll","NotAll")) %>% 
  mutate(perc = value/92 * 100) %>% 
  ggplot(., aes(x = dataset, y = perc)) + geom_bar(stat = "identity") +
  mytheme + labs(x = "", y = "ERCCs detected (%)")


p1 

Redundant Isoforms from Chained Data

Plot is drawn from All_Merged dataset in Sqanti chained output

redundant <- unique[,c("FL.All_Merged")] %>% table() %>% melt 
colnames(redundant) <- c("ERCC","num_isoforms")

p2 <- redundant %>% 
  group_by(num_isoforms) %>% tally() %>% ggplot(., aes(x = num_isoforms, y = n)) + geom_bar(stat = "identity") + 
  mytheme + labs(x = "ERCC", y = "Number of Isoforms")

p2

datatable(redundant)

Isoforms vs Concentration

isoform_conc <- merge(redundant, ERCC_conc, by.x = "ERCC", by.y = "ERCC_ID", all = TRUE) %>%
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  replace_na(list(num_isoforms = 0)) %>% 
  mutate(num_isoforms = as.factor(num_isoforms))

p3 <- ggplot(isoform_conc, aes(x = num_isoforms, y = log2_amount_of_ERCC, colour = num_isoforms)) + geom_jitter(width = 0.2) + theme(legend.position = "none") + 
  mytheme + labs(x = "Number of Isoforms", y = "Amount of ERCC (log2)")

p3

Pipeline

Pipeline for analysis

Classification file

# Read in sqanti classification file mapped to ERCC (post chaining WT and TG)
class <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/SQANTIAll/All_Merged.collapsed.filtered_classification.filtered_lite_classification.txt", header = T, as.is = T, sep = "\t")

datatable(head(class), options = list(scrollX = TRUE))

Unique ERCC

cat("Total unique ERCCs:", length(unique(class$chrom)), paste0("(", round(length(unique(class$chrom))/92 * 100,2), "%)"))

## Total unique ERCCs: 37 (40.22%)

Redundant Isoforms

redundant <- class[,c("chrom")] %>% table() %>% melt 
colnames(redundant) <- c("ERCC","num_isoforms")

p4 <- redundant %>% 
  group_by(num_isoforms) %>% tally() %>% ggplot(., aes(x = num_isoforms, y = n)) + geom_bar(stat = "identity") + 
  mytheme + labs(x = "ERCC", y = "Number of Isoforms")

p4

datatable(redundant)

Isoforms vs Concentration

isoform_conc <- merge(redundant, ERCC_conc, by.x = "ERCC", by.y = "ERCC_ID", all = TRUE) %>%
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  replace_na(list(num_isoforms = 0)) %>% 
  mutate(num_isoforms = as.factor(num_isoforms))

p5 <- ggplot(isoform_conc, aes(x = num_isoforms, y = log2_amount_of_ERCC, colour = num_isoforms)) + geom_jitter(width = 0.2) + theme(legend.position = "none") + 
  mytheme + labs(x = "Number of Isoforms", y = "Amount of ERCC (log2)")

p5

Pipeline

Pipeline for analysis

Classification file

class <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/SQANTIDEMUX/All_Merged.collapsed.filtered_classification.filtered_lite_classification.txt",header = T, as.is = T, sep = "\t") %>% 
  mutate(Total_FL = .[,46:55] %>% rowSums(na.rm = TRUE))
 
datatable(head(class), options = list(scrollX = TRUE))

Unique ERCC

cat("Total unique ERCCs in both WT and TG:", length(unique(class$chrom)), paste0("(", round(length(unique(class$chrom))/92 * 100,2), "%)"))

## Total unique ERCCs in both WT and TG: 37 (40.22%)

# replace FL count with pb id, otherwise NA if 0
samples <- colnames(class)[46:55]
unique <- setNames(data.frame(matrix(ncol = 12, nrow = nrow(class))), c("isoform","ERCC",samples))
unique$ERCC <- class$chrom
unique$isoform <- class$isoform
for (j in samples){
  for(i in 1:nrow(class)){unique[[j]][i] <- ifelse(class[[j]][i] > 0, paste(class$chrom[i]),paste("NA"))}
}

# number of ERCC detected across WT and TG 
totalERCC <- apply(unique[,-c(1,2)], 2, function(x) length(unique(x))) %>% melt() %>% 
  rownames_to_column(., "dataset") %>% 
  mutate(phenotype = word(dataset, c(2), sep = "_")) %>%
  `colnames<-`(c("sample", "num", "phenotype"))

aggregate(x = totalERCC$num, by = list(totalERCC$phenotype),FUN = mean) %>% `colnames<-`(c("Phenotype", "Mean_ERCC_Detected")) %>%
  mutate(perc = round(Mean_ERCC_Detected/92 * 100,2))

##   Phenotype Mean_ERCC_Detected  perc
## 1        TG               32.2 35.00
## 2        WT               32.4 35.22

# count the unique pbid (/ERCC) per dataset
p6 <- apply(unique[,-c(1,2)], 2, function(x) length(unique(x))) %>% melt() %>% 
  rownames_to_column(., "dataset") %>%
  mutate(Genotype = word(.$dataset,c(2), sep = "_")) %>% 
  mutate(perc = value/92 * 100) %>% 
  ggplot(., aes(x = Genotype, y = perc)) + geom_boxplot() +
  geom_jitter(aes(color = Genotype), shape=16, position=position_jitter(0.2)) + 
  mytheme + labs(x = "", y = "ERCCs detected (%)", title = "\n\n") +
  theme(legend.position = "none") + ylim(20,60) +
  scale_color_manual(values = c(label_colour("TG"),label_colour("WT")))


p6

Redundant Isoforms

redundant <- unique[,c("ERCC")] %>% table() %>% melt 
colnames(redundant) <- c("ERCC","num_isoforms")

p7 <- redundant %>% 
  group_by(num_isoforms) %>% tally() %>% ggplot(., aes(x = num_isoforms, y = n)) + geom_bar(stat = "identity") + 
  mytheme + labs(x = "ERCC", y = "Number of Isoforms")

p7

datatable(redundant)

# number of ERCC with redudant isoforms
cat("Number of ERCC with more than one isoform", nrow(redundant[redundant$num_isoforms > 1,]), 
    "(",round(nrow(redundant[redundant$num_isoforms > 1,])/92 * 100, 2),"%)")

## Number of ERCC with more than one isoform 8 ( 8.7 %)

Isoforms vs Concentration

isoform_conc <- merge(redundant, ERCC_conc, by.x = "ERCC", by.y = "ERCC_ID", all = TRUE) %>%
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  replace_na(list(num_isoforms = 0)) %>% 
  mutate(num_isoforms = as.factor(num_isoforms))

p8 <- ggplot(isoform_conc, aes(x = num_isoforms, y = log2_amount_of_ERCC, colour = num_isoforms)) + geom_jitter(width = 0.2) + theme(legend.position = "none") + 
  mytheme + labs(x = "Number of Isoforms", y = "Known Amounts of ERCC (log2)", title = "\n\n")

p8

+ TAMA

To use TAMA’s tama_remove_fragment_models.py to remove shorter redundant “isoforms” from the same ERCCs. Script is applied after SQANTI filtering.

TAMA <- "/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/TOFU/TAMA_Filter/"
tama_retained <- read.table(paste0(TAMA,"All_Merged_postsqanti.bed")) %>% mutate(TAMA = "Retained") 
tama_discard <- read.table(paste0(TAMA,"All_Merged_postsqanti_discarded.txt")) %>% mutate(TAMA = "Discarded")
tama <- bind_rows(tama_retained, tama_discard) %>% mutate(isoform = word(.$V4, c(2),sep = ";"))

p9 <- tama %>% group_by(V1,TAMA) %>% tally() %>%
  ggplot(., aes(x = reorder(V1,-n), y = n, fill = TAMA)) + geom_bar(stat = "identity") + 
  labs(y = "Number of Isoforms", x = "", title = "\n\n") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.title = element_blank(), legend.position = c(0.9, 0.9)) +
  mytheme 

# reads filtered from TAMA 
# total FL reads across 10 samples associated with isoform 
p10 <- class[class$chrom %in% tama_discard$V1,] %>%
  mutate(Total_FL = .[, 46:55] %>% rowSums(na.rm = TRUE)) %>%
  left_join(., tama, by = ("isoform")) %>% 
  mutate(log10_Total_FL = log10(Total_FL)) %>%
  ggplot(., aes(x = chrom, y = log10_Total_FL, colour = TAMA)) + geom_point() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + labs(y = "Full-Length Reads (Log10)", x = "", title = "\n\n") 

p11 <- class %>%
  mutate(Total_FL = .[, 46:55] %>% rowSums(na.rm = TRUE)) %>%
  left_join(., tama, by = ("isoform")) %>% 
  mutate(log10_Total_FL = log10(Total_FL)) %>%
  ggplot(., aes(x = chrom, y = log10_Total_FL, colour = TAMA)) + geom_point() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + labs(y = "Full-Length Reads (Log10)", x = "") 

p9 

p10

Correlation

ERCC_corr <- tama_retained %>% 
  mutate(isoform = word(V4, c(2),sep = ";")) %>%
  left_join(.,class[,c("isoform","Total_FL")], by = "isoform") %>% 
  left_join(.,ERCC_conc, by = c("V1" = "ERCC_ID")) %>% 
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  mutate(log2_FL_reads = log2(Total_FL)) 

p12 <- density_plot(ERCC_corr,"log2_amount_of_ERCC","log2_FL_reads", "Amount of ERCC (Log2)", "Number of Full-Length Reads (Log2)","")

## 
##  Pearson's product-moment correlation
## 
## data:  dat[[x.var]] and dat[[y.var]]
## t = 8.5057, df = 35, p-value = 4.89e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6770173 0.9043953
## sample estimates:
##       cor 
## 0.8209479 
## 
## [1] ""
## [1] "Correlation betweenlog2_amount_of_ERCCandlog2_FL_reads"
## [1] "corr.value0.820947870699425"
## [1] "p.value4.89015691673981e-10"

p12

Alignment vs Identity of all mapped ERCCs

Threshold:
minimum alignment coverage = 0.95 (parameter set from discarded ERCC)
minimum alignment identity = 0.95

cat("Total unique ERCCs in both WT and TG:", length(unique(Merged_mapping$ERCC_mapped$V2)), paste0("(", round(length(unique(Merged_mapping$ERCC_mapped$V2))/92 * 100,2), "%)"))

## Total unique ERCCs in both WT and TG: 58 (63.04%)

# plot of Alignable length and identity
p13 <- Merged_mapping$ERCC_mapped %>% 
  mutate(identity = V8 * 100) %>% 
  mutate(length = V6 * 100) %>%
  ggplot(., aes(x = length, y = identity)) + 
  geom_point() + mytheme + labs(x = "Alignment Length (%)", y = "Alignment Identity (%)", title = "\n\n") +
  geom_hline(yintercept = 95, linetype = 2) +
  geom_vline(xintercept = 95, linetype = 2) +
  annotate("text", x = 20, y = 94, label = "95% Threshold") +
  annotate("text", x = 94, y = 50, label = "95% Threshold", angle = 90) +
  annotate("rect", xmin = 95, xmax = 100, ymin = 95, ymax = 100, fill = "palegreen", alpha = 0.2)

p13 

Post Filtering

filter_map <- Merged_mapping$ERCC_mapped[Merged_mapping$ERCC_mapped$V8 > 0.95 & Merged_mapping$ERCC_mapped$V6 > 0.90,]
cat("Total unique ERCCs in both WT and TG (post filtering on length and identity):", length(unique(filter_map$V2)), paste0("(", round(length(unique(filter_map$V2))/92 * 100,2), "%)"))

## Total unique ERCCs in both WT and TG (post filtering on length and identity): 57 (61.96%)

cluster <- filter_map %>% group_by(V2) %>% tally() %>% 
  mutate(V1 = V2) %>%
  full_join(., tama[tama$TAMA == "Retain",c("V1","TAMA","isoform")], by = "V1") %>% 
  .[,c("V2","n","TAMA","isoform")] %>%
  `colnames<-`(c("ERCC_ID", "Num_Isoform", "Status","isoform")) %>% 
  full_join(.,ERCC_conc[,c("ERCC_ID","amount_of_ERCC")], by = "ERCC_ID") %>%
  replace_na(list(Num_Isoform = 0, Status = "Discard")) %>%
  mutate(Status=replace(Status, Num_Isoform==0, "Not_Detected"))

p14 <- ggplot(cluster,aes(x = Num_Isoform, y = log2(amount_of_ERCC), colour = Status)) + geom_jitter(width = 0.2) +
  mytheme + labs(x = "Number of Isoforms", y = "Amount of ERCC (log2)") +
  scale_colour_discrete(name = "", label = c("Discard","Not Detected", "Retain")) + 
  theme(legend.position = c(0.8, 0.2))

p14

Discarded ERCC

Rationale for discarded ERCC (1 isoform) from cupcake collapse
Cupcake’s parameters: minimum alignment coverage = 0.99, minimum alignment identity = 0.95

# rationale for discarded ERCC (1 isoform)
discarded_ERCC <- cluster[cluster$Status == "Discard" & cluster$Num_Isoform == "1",]
discarded_transcript <- filter_map[filter_map$V2 %in% discarded_ERCC$ERCC_ID,]

cupcake_discard <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/TOFU/All_Merged.ignored_ids.txt", as.is = T, sep = "\t") %>% `colnames<-`(c("transcript", "status"))

p15 <- cupcake_discard %>% filter(transcript %in% discarded_transcript$V1) %>% 
  left_join(., discarded_transcript[,c("V1","V2")], by = c("transcript" = "V1")) %>% 
  mutate(coverage = as.numeric(word(status, c(2), sep = " ")) * 100) %>%
  ggplot(., aes(x = V2, y = coverage)) + geom_point() + ylim(80,100) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + labs(x = "Discarded ERCC", y = "Alignment Length (%)")

p15

datatable(cupcake_discard %>% filter(transcript %in% discarded_transcript$V1))

Unique ERCC

Reduced cupcake’s collapse default parameters (0.99) to 0.95

class <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/TOFU_ADJUST/SQANTI/All_Merged.collapsed.filtered_classification.filtered_lite_classification.txt",header = T, as.is = T, sep = "\t") 

class <- class %>% mutate(Total_FL = .[,46:55] %>% rowSums(na.rm = TRUE))

cat("Total unique ERCCs in both WT and TG:", length(unique(class$chrom)), paste0("(", round(length(unique(class$chrom))/92 * 100,2), "%)"))

## Total unique ERCCs in both WT and TG: 57 (61.96%)

# replace FL count with pb id, otherwise NA if 0
samples <- colnames(class)[46:55]
unique <- setNames(data.frame(matrix(ncol = 12, nrow = nrow(class))), c("isoform","ERCC",samples))
unique$ERCC <- class$chrom
unique$isoform <- class$isoform
for (j in samples){
  for(i in 1:nrow(class)){unique[[j]][i] <- ifelse(class[[j]][i] > 0, paste(class$chrom[i]),paste("NA"))}
}


# count the unique pbid (/ERCC) per dataset
p16 <- apply(unique[,-c(1,2)], 2, function(x) length(unique(x))) %>% melt() %>% 
  rownames_to_column(., "dataset") %>%
  mutate(Genotype = word(.$dataset,c(2), sep = "_")) %>% 
  mutate(perc = value/92 * 100) %>% 
  ggplot(., aes(x = Genotype, y = perc)) + geom_boxplot() +
  geom_jitter(aes(color = Genotype), shape=16, position=position_jitter(0.2)) + 
  mytheme + labs(x = "", y = "ERCCs detected (%)", title = "\n\n") +
  theme(legend.position = "none") + 
  ylim(20,60) +
  scale_color_manual(values = c(label_colour("TG"),label_colour("WT")))


p16

Redundant Isoforms

redundant <- unique[,c("ERCC")] %>% table() %>% melt 
colnames(redundant) <- c("ERCC","num_isoforms")

p17 <- redundant %>% 
  group_by(num_isoforms) %>% tally() %>% ggplot(., aes(x = num_isoforms, y = n)) + geom_bar(stat = "identity") + 
  mytheme + labs(x = "ERCC", y = "Number of Isoforms")

p17

datatable(redundant)

Isoforms vs Concentration

isoform_conc <- merge(redundant, ERCC_conc, by.x = "ERCC", by.y = "ERCC_ID", all = TRUE) %>%
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  replace_na(list(num_isoforms = 0)) %>% 
  mutate(num_isoforms = as.factor(num_isoforms))

p18 <- ggplot(isoform_conc, aes(x = num_isoforms, y = log2_amount_of_ERCC, colour = num_isoforms)) + geom_jitter(width = 0.2) + theme(legend.position = "none") + 
  mytheme + labs(x = "Number of Isoforms", y = "Amount of ERCC (log2)")

p18

+ TAMA

To use TAMA’s tama_remove_fragment_models.py to remove shorter redundant “isoforms” from the same ERCCs. Script is applied after SQANTI filtering.

TAMA <- "/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/IsoSeq/Whole_Transcriptome/All_Samples/ERCC/TOFU_ADJUST/TAMA_Filter/"
tama_retained_adjusted <- read.table(paste0(TAMA,"All_Merged_postsqanti.bed")) %>% mutate(TAMA = "Retain") 
tama_discard <- read.table(paste0(TAMA,"All_Merged_postsqanti_discarded.txt")) %>% mutate(TAMA = "Discard")
tama <- bind_rows(tama_retained_adjusted, tama_discard) %>% mutate(isoform = word(.$V4, c(2),sep = ";"))

p19 <- tama %>% group_by(V1,TAMA) %>% tally() %>%
  ggplot(., aes(x = reorder(V1,-n), y = n, fill = TAMA)) + geom_bar(stat = "identity") + labs(y = "Number of Isoforms", x = "ERCC") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme

# reads filtered from TAMA 
# total FL reads across 10 samples associated with isoform 
p20 <- class[class$chrom %in% tama_discard$V1,] %>%
  mutate(Total_FL = .[, 46:55] %>% rowSums(na.rm = TRUE)) %>%
  left_join(., tama, by = ("isoform")) %>% 
  mutate(log10_Total_FL = log10(Total_FL)) %>%
  ggplot(., aes(x = chrom, y = log10_Total_FL, colour = TAMA)) + geom_point() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + labs(y = "FL Reads (Log10)", x = "ERCC") 

p11 <- class %>%
  mutate(Total_FL = .[, 46:55] %>% rowSums(na.rm = TRUE)) %>%
  left_join(., tama, by = ("isoform")) %>% 
  mutate(log10_Total_FL = log10(Total_FL)) %>%
  ggplot(., aes(x = chrom, y = log10_Total_FL, colour = TAMA)) + geom_point() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + labs(y = "FL Reads (Log10)", x = "ERCC") 

p19 

p20

Rescued lowly-expressed transcripts

#setdiff(tama_retained$V1,tama_retained_adjusted$V1)
#setdiff(tama_retained_adjusted$V1,tama_retained$V1)

tama_retained_adjusted$status <- ifelse(tama_retained_adjusted$V1 %in% tama_retained$V1, "99% Identity", "95% Identity")

p21 <- tama_retained_adjusted %>% 
  left_join(.,ERCC_conc, by = c("V1" = "ERCC_ID")) %>%
  ggplot(., aes(x = status, y = log2(amount_of_ERCC))) + geom_boxplot() + geom_jitter(width = 0.2) + 
  mytheme + labs(x = "", y = "Known Amounts of ERCC (log2)", title = "\n\n")

p21

Correlation

ERCC_corr <- tama_retained_adjusted %>% 
  mutate(isoform = word(V4, c(2),sep = ";")) %>%
  left_join(.,class[,c("isoform","Total_FL")], by = "isoform") %>% 
  left_join(.,ERCC_conc, by = c("V1" = "ERCC_ID")) %>% 
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  mutate(log2_FL_reads = log2(Total_FL)) 

p22 <- density_plot(ERCC_corr,"log2_amount_of_ERCC","log2_FL_reads", "Amount of ERCC (Log2)", "Number of Full-Length Reads (Log2)","")

## 
##  Pearson's product-moment correlation
## 
## data:  dat[[x.var]] and dat[[y.var]]
## t = 38.688, df = 55, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9697060 0.9894717
## sample estimates:
##      cor 
## 0.982118 
## 
## [1] ""
## [1] "Correlation betweenlog2_amount_of_ERCCandlog2_FL_reads"
## [1] "corr.value0.982117953981347"
## [1] "p.value1.41213399302975e-41"

p22

# outputs
pdf(paste0(output_plot_dir,"/ERCC.pdf"), width = 15, height = 11.5)
plot_grid(p6,p16,p12,p22, labels = c("a","b","c"), label_size = 25, label_fontfamily = "CM Roman", label_fontface = "bold",scale = 0.95)
top_row <- plot_grid(p8,p10, labels = c("a","b"), label_size = 25, label_fontfamily = "CM Roman", label_fontface = "bold")
plot_grid(top_row,p9, labels = c('', 'c'), label_size = 25, ncol = 1,label_fontfamily = "CM Roman", label_fontface = "bold", scale = 0.95)
plot_grid(p13,p21, labels = c("a","b"), label_size = 25, label_fontfamily = "CM Roman", label_fontface = "bold",scale = 0.95, ncol = 1)
plot_grid(pmousemap,pERCCmap, labels = c("a","b"), label_size = 25, label_fontfamily = "CM Roman", label_fontface = "bold",scale = 0.95)
dev.off()

## png 
##   2